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We extend the computational visualization method proposed in PJ |2] using the 
concept of the harmonic time averages presented in j3]. Algorithms for frequency 
analysis of the phase space are constructed and implemented numerically, pro- 
ducing a graphical visualization of the periodic sets for a given periodicity. The 
convergence of the harmonic time averages is addressed, and an algorithm using 
more functions is proposed for visualization of the phase space periodic partitions. 
Visualization of chaotic regions based on the same concept is exposed as well. The 
method is presented in the context of the discrete-time dynamical systems using the 
Chirikov standard map as the prototype. Applications to other maps are included. 
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1 Introduction 

Methods of computational investigation of complex dynamical systems are of crucial importance today, 
given the increasing range of interdisciplinary phenomena currently under examination [4 a . As virtually 
none of the naturally motivated dynamical system can by entirely treated with analytical techniques, the 
necessity of numerical investigation is beyond doubt. When choosing from a diverse spectrum of available 
computational approaches, one takes into account the nature of the problem as well as the study direction 
that is to be followed. 

Typically, when the dynamical systems are considered, the direct integration of system's trajectories 
is performed for long integration times, seeking to optimize the precision with respect to the numerical 
cost. The obtained trajectories are then analyzed, either in the context of specific trajectory-approach 
to the system, or by looking for their universal properties that report about global nature of the system 
under investigation. In particular, one is often interested in frequency- analysis of the trajectory done by 
computing the power spectrum of the trajectory's time-signal. Power (Fourier) spectrum of a time-signal 
gives a frequency decomposition of the trajectory, showing it in terms of its frequencies [5j |6] . 

The idea of frequency analysis of time series inspired new approaches to study of the dynamical systems 
[7] , specifically in the context of chaotic motion analysis [HI IH] > high-dimensional dynamical systems [TU] 
and the study of celestial dynamics 11J. Namely, an algorithm can be constructed based on generalization 
of power spectrum analysis that estimates the system's fundamental frequencies in a decreasing order [5] . 
Approaches of this sort were very useful, both in the context of continuous-time [7] and discrete-time [12] 
dynamical systems. More recently, a new approach has been proposed based on wavelet frequency analysis 



that investigates the properties of specific chaotic orbits and detects resonance trappings and transitions 
[13U14j . In particular, the method is suited for examination of weakly chaotic orbits. Furthermore, another 
recent time-frequency analysis method for chaotic time series based on the power spectrum estimator was 
suggested, that distinguishes between noisy flow and the colored noise [15] . 

The key downside of the mentioned approaches lies in their locality: the frequency analysis is usually 
performed on single trajectories/time-signals only. While providing a detailed analysis of trajectories 
under examination, the method remains focused on pre-defined ensembles of trajectories, unable to give 
conclusions regarding the global characteristics of the system, that would follow from comparison among 
frequency properties of many trajectories. In the present work, we extend the previously exposed idea of 
invariant sets visualization method based on ergodic partitioning [I], to the domain of global frequency 
analysis of the dynamical phase space. We construct harmonic time averages of functions defined over 
the phase space, and produce color-plots reporting about frequency phase space structure by visualizing 
periodic subsets of prescribed frequencies. We therefore obtain algorithm that produces detailed analysis 
of phase space in terms of periodic sets resonating with a single frequency: by using more frequencies we 
extend the method to full phase space frequency analysis. By employing parallel computing we are able to 
save numerical time and obtain global phase space plots relatively fast. 

As opposed to other approaches, our method is global in its nature, as it simultaneously provides insights 
into the entire phase space. Also, along the lines of our previous study [T], this method is suitable for high- 
dimensional systems in terms of considering phase space sections of smaller dimensionality. Furthermore, 
as it will be shown, the method is applicable to visualization of the chaotic zone of the dynamical system. 
We will examine the method using well-known example of standard map, as this systems possesses a large 
variety of periodic sets and a chaotic zone whose locations are known. 

This paper is organized as follows: we briefly expose the theoretical background of our method in Section 
[2] and show its implementation with single function in Section [3] The convergence and precision issues 
are settled in Section [4] In Section [5] we show the application of multi-functional approach, designing a 
simple algorithm for visualization of periodic and chaotic partitions. We show the applicational extent of 
out method in Section [6] , and give a few concluding remarks in Section [7] 



2 The Visualization Idea 

In this Section we present the concept of harmonic time averages [5] that relate the ergodic theory ideas 
to the frequency analysis. We expose the basics of the mathematical background necessary for our study; 
for more rigorous details we however address the reader to [3j I16[ 117] . 

Let a discrete-time measure-preserving map x' = Tx on a compact phase space A C K™ be denoted as: 

x n+ i = Tx„ or x„ = T"x n > 0. (1) 

Our central aim is to visualize the periodic sets B p C A with the periodicity p for the map T, defined as: 

x„eB p => T pxn x eB p Vn>0. (2) 

Periodic sets are essentially generalizations of the periodic orbits, with the period-1 set being the usual 
invariant set. Note that there are always p distinct period-p sets B^,.. such that a phase space point x 
makes p jumps among them before completing a whole cycle. Union of a sequence of such sets UfeBp = B 
is an invariant set composed of p disjoint parts, that we will call a periodic chain. In what follows we will 
expose a computational method based on ergodic theory designed to graphically visualize such set in the 
dynamics phase space, based on the methods described in [T]. 

Consider L 2 real-valued functions on A and let the harmonic time average /*(xq) of a function / G 
L 2 (A) for the frequency u) G [0, |] corresponding to a phase space point xo G A, be defined as: 

n— 1 

/:(x)= lim -5> l2 ^/(T fc x). (3) 

fc=0 



2 



By the Ergodic Theorem this limit exists almost everywhere in the phase space for any measure-preserving 
map T |161 1171 118j . Harmonic time averages are a generalization of the time averages known from the 
standard ergodic theory , which is easy to see by observing that /* =0 = /* |19j . Furthermore, from the 
definition Eq. ([3| it follows: 

£(T»x) = e - 42 ™"/:(x), (4) 

which hence implies: 

|£(T»x)| = |/:(x)| Vn>0, (5) 

meaning that while f* is not an invariant function, its absolute value \f*\ is. That is to say each trajectory 
(periodic or not) has a constant |/^| for each point. As throughout this work we will be considering only 
the absolute values of the harmonic time averages, we introduce the notation: 

We will refer to functions /i w (x) as harmonic time averages, despite really intending their absolute values. 
Consider now the circle S 1 = [0, 2ir[ and the shift map W mapping the circle onto itself: 

0^6 = 6 + 2™, (6) 

for some constant angle u> € [0, |]. Given a map T : A — » A and a set B c A, a shift map Q u is called the 
factor map to T on B if there is a measure-preserving homeomorphism F : A — > S 1 such that: 

(FoT)(x) = (e B oF)(x), (7) 

Vx £ B. This means the dynamics T on B C A is topologically equivalent to a shift map, which implies 
that the set B is a periodic set or a union of periodic sets with the period 1/lo. The following results hold: 



Theorem Given a dynamics T : A — > A and / £ L 2 (A), if the harmonic time average h u is non-zero on 
some set B C A then there exists a factor map for T on the set B given by the shift map with the angle lo. 



Theorem If T : A — ► A admits a shift map with the angle lo as a factor map on some set B C A, then 
there exists an / S L 2 (A) such that its harmonic time average 7i w is non-zero on _B. 

Therefore, the harmonic time averages for the specific frequencies u> can be used for detecting the periodic 
sets and periodic chains of the period l/u>. If h^'s frequency ui resonates with the frequency of a periodic 
set, its final limit value for the points within this set will be non-zero; otherwise the summation of the 
complex phases e _l27r ' im will average out its limit value to zero. Computing the harmonic time average 
over the whole phase space will therefore expose only the phase space subsets whose frequency resonates 
with lo, otherwise h u will average out to zero (we will define the periodic sets mostly by their frequency 
that is inverse of their periodicity) . 

Observe the relationship between a harmonic time average and a Fourier transform: while the later 
gives the entire frequency spectrum of a certain time-signal (a single trajectory), the former gives the 
phase space partition that resonates with the particular chosen frequency. A harmonic time average value 
for a given trajectory is simply the Fourier transform's value at that frequency for this trajectory. While 
the Fourier transform provides us with a full frequency analysis for a single trajectory, the harmonic time 
average gives a global phase space analysis, but for a single frequency only. 

Furthermore, a shift map with an angle lo is also a shift map with all the angles that are integer multiples 
of lo. Therefore, a harmonic time average of frequency lo reveals as non-zero all the periodic sets with the 
periodicities that are integer multiples of 1/lo: a harmonic time average for lo — \ will not only resonate 
with all the period-2 set, but also all the even-period sets. This reduced the sequence of interesting periods 
to the numbers that are mutually prime. Throughout this work we will consider periodicities given by the 
sequence of prime numbers: 2,3,5,7 etc. We will not consider the period-1 sets (frequency lo = 0) as /i w= o 
will simultaneously expose periodic sets of all other periodicities (as expected given that h u = = |/*|), 
which does not allow a detailed analysis. 
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The final value of a harmonic time average over a resonating periodic set will also be determined by the 
properties of the function / in question: the limit value for a resonating periodic set can happened to be 
zero due to location of the set itself with respect to the properties of the function. For this reason we will, 
in a way similar to our previous study PQ, employ the harmonic time averages of more functions for the 
same frequency in order to optimize the visualization of periodic sets of that frequency. The computational 
algorithm and its implementation are described in the rest of this Section. 

2.1 The computing algorithm and the numerical details 

We construct the algorithm for the phase space partitioning into a periodic partition of a given frequency 
using harmonic time averages of multiple functions. For example, the periodic partition for the frequency 
co = | is constructed as a partition of all even-period periodic sets (and of course the non-resonating rest 
of the phase space) . Following [T] , we limit the discussion to the case of a two-dimensional map with the 
phase space A — [0, 1] x [0, 1]. 

step 1 Set up a grid (lattice) of initial grid-points (xo, yo) on the phase space A and select the visualization 
frequency uj in form ^ with p prime number 2,3,5,7 etc. 

step 2 Pick N linearly independent functions {/i, . . . Jn} from L 2 (A) and calculate their partial harmonic 
time averages for ng na j iterations for each initial grid-point, which approximate the real harmonic 
time averages {h^, . . . } 

step 3 To every initial grid-point (xq, yo) associate the harmonic time average vector corresponding to it: 
(xo,Uo) — ► ^uj(x Q ,y ) = {hl(x ,yo),...,h^(xo,y Q )} S [0, 1} N (8) 

step 4 Observe the distribution of the harmonic time average vectors h w (xo, yo) throughout M. N and group 
them optimally into clusters. Divide A into a union of subsets, with each subset being given by those 
grid-points (xo,yo) whose harmonic time average vectors h u (xo,yo) belong to the same cluster of 
vectors. This union of subset is the iV-order approximation periodic partition of A with frequency lu. 

The optimal number of iterations ng n „j is to be set in accordance with the harmonic time averages con- 
vergence properties that will be studied in the Section [4] Visualization also depends on the clustering 
procedure, that we will examine in the Section [5] The choice of functions {/i, . . . /jv} is again done by 
looking at an orthogonal basis on L 2 (A) as it provides the simplest source of linearly independent functions. 
Throughout this work we will employ functions from the 2D Fourier orthogonal basis given as: 

sin(7rmc + irmy) with n, m G Z. (9) 

These functions are among them linearly independent for any n =/= m. 

In a way identical to pQ we will be using parallel processing to enhance the efficiency of computation. 
We typically use the grid of 800 x 800 initial grid-points and iterate the dynamics for ng na ^ = 30,000 
iterations for each grid-point to approximate the harmonic time average value for each function. It takes 
about 10-15 minutes to compute one time average on a single processor for a grid of this size, a typical run 
was however done using 5-10 processors (depending on the availability) and took about 5 minutes for a 
single function. Computing more functions simultaneously insignificantly increases the computation time. 

2.2 The standard map as the testing prototype 

As in the previous work, we will rely on the Chirikov's standard map |20j for testing the described method, 
as this map possesses a variety of periodic chains and sets and other dynamical behaviors for a wide range 
of parameters [21. 5 22J. The map is a homeomorphism on a 2D torus given by: 

x 1 = x + y + £sin(27ra) [modi] 
y' = y + esin(27rx) [modi] 

where (x,y) £ [0, 1] x [0, 1] = [0, l] 2 (the usual standard map's parameter k is here k — 2ire). 
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3 Phase Space Plots with a Single Function 



In this Section we consider the simplest case of N = 1 and show the color-plots of harmonic time averages 
for a single function under the dynamics of the standard map Eq. (10 1. As mentioned, we are using a grid 



of 800 x 800 initial grid-points and n.g na j = 30, 000 iterations. For simplicity reasons, in this Section we 
will focus only on the perturbation value e — 0.12 and the function f(x,y) — sin(27rizi + 3iry). 

In Figs.[l]i,b&c we show plots for the frequencies uj — |, |, | that are visualizing the periodic sets with 
the periodicities that are multiples of 2, 3 and 5 respectively. For clarity we are using a log-scale in for the 
hu value and stretching the colorbar from the smallest to the biggest value of log h w . The chains (families) 





(a) 



(b) 




| 



(c) 

Figure 1: Harmonic time averages \ogh u of sin(27ra; + iity) for the map 
frequencies ui = ^ in (a), u> = j in (b), ui = | in (c) and uj = — in (d). 



(d) 

Eq. (fill with £ = 0.12 with the 



of periodic sets are correctly predicted in all the plots, together with the periodic sets of higher (integer 
multiple) periodicities. The final harmonic time average value h u is modulated by the values that the 
employed function is taking over the the periodic set in question, which explains why are some periodic 
chains better visible than the others. Recall that the large period-2 island around the elliptic fixed point 
(§,5) (cf. Fig.jT^) is actually a nested set of infinitely many quasi-periodic orbits enclosing the fixed point. 
Each of these quasi-periodic orbits is of course a period-2 set, but the color-difference among them is less 
visible due to the use of the log-scale. By using the log-scale we are focusing on revealing the periodic 
chains themselves as phase space subsets, rather then trying to discern their internal sub-structure. 

Some of the chains of periodic sets are visible in more plots (like period-6 chain in Fig.[l^L and Fig. [I]}), 
as their periodicities are common-multiples of various basic frequencies. On the other hand, some periodic 
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chains (e.g. period-10 chain) are visible in one plot but not in the others (Fig.[Tp and Fig.[Tp) despite 
resonating with both frequencies: as already mentioned, this is due to the properties of the function used 
for averaging, and will be taken care of by using multiple functions (see the forthcoming Sections). 

Also, observe that the chaotic zone (already locally present for this e-value around the hyperbolic fixed 
point (0,0) = (1,1)) is resonating at all frequencies, although much more weakly. Moreover, the chaotic 
zone is also weakly resonating in case of an irrational frequency 1/n as visible in Fig.[T]i, at which no other 
periodic set resonates at all (within the limits of the precision of numerically irrational number). As we 
shall see, this can be used for visualization of the chaotic zone itself using more functions. 



4 The Convergence Study 

In this Section we study the convergence of the harmonic time averages addressing the questions of the 
precision of their final values and the optimal number of iterations required. We will distinguish between 
the regular (periodic) and chaotic orbits, and between rational and irrational frequencies, justifying the 
visualization results from the previous Section. 

Consider the n-th partial harmonic time average for some frequency u> given by: 

n— 1 

K(xo,Vo) = -£e"/(T fc (z 0)2/0 )), (11) 

fe=0 

with lhrin^oo |/"(^Oj2/o)| — ^(^cbyo) assumed to exist for all the grid-points (xo,yo)- The difference: 

A (x , V0 )(n) = | K(x ,y ) - \fZ(x ,y )\ | (12) 

in function of n is a sequence whose behavior is to be studied in relation to the orbits/frequencies mentioned 
above. We iterate the dynamics for 10 7 iterations recording the first 10 6 iterations and defining | = 

h u . This gives the sequence A(n) whose asymptotic properties are investigated. 

The Regular Region. As mentioned earlier, for all the regular points (trajectories), the harmonic time 
average converges to zero in case its frequency mismatches the periodicity of the underlying set, and 
otherwise converges to some positive value smaller then 1. In analogy to what observed in p~j[23], the 
convergence pattern of A(n) in this region neatly follows 1/n regardless of the frequency, as shown in the 
convergence plot Fig. [2^. Moreover, this holds also in the case of irrational frequency, as illustrated in 
Fig.[2|) (all h^s converge to zero). These results apply universally for all the regular points/periodic sets 
independently from the e- value, periodicity or the function involved, allowing a good precision estimation. 

The Chaotic Region. In further analogy to the previous studies 1 , 23 , the convergence pattern exhibited 
by A(n) in this region is at best given by l/-\/ n an d hence bounded by a — | slope, which is in general better 
pronounced at bigger e- values. At smaller e- values A(n) is at best asymptotically bounded by the weak 
slopes between —1 and zero (see discussion of weakly chaotic region in [T]). These results hold uniformly 
for both rational and irrational frequencies, as shown in Figs.[2j;&d regardless of the particular grid-point, 
e-value or function involved. 

Note that in analogy with the case of the usual time averages, the precision of (weakly) chaotic trajec- 
tories is to determine the lower precision bound of a harmonic time average phase space plot. We therefore 
set the iteration limit to «g na i = 30, 000 having the overall precision of O(10 -2 ) or better. 

Furthermore, observe that in the case of an irrational frequency, a harmonic time average always 
converges to zero, but with a rate that depends on the nature of the underlying trajectory. For a periodic 
set of any periodicity, this limit will converge to zero in way identical to the case of a harmonic time average 
with a non-matching frequency, whereas for any chaotic orbit it will always "weakly resonate" converging 
at a slower rate of A(n) ~ This provides a simple way to approximately visualize the chaotic zone as 
it will be shown later in this work. 
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5 The Phase Space Plots with Multiple Functions 

In this Section we consider the harmonic time averages corresponding to more linearly independent func- 
tions for the same frequency. In analogy with the construction of the ergodic partition pQ , we examine the 
correspondence given by the Eq. ^ and seek conclusions regarding the periodic sets phase space distribu- 
tion by investigating clustering of the harmonic time average vectors. We also propose a simple algorithm 
for visualizing the periodic partition for a given frequency to and the phase space ergodic regions. 

5.1 Two-function scatter plots 

We start by TV = 2 case considering two linearly independent functions for a fixed frequency ur. 

(a;o,yo)egrid — > {hi, h 2 J G [0, l] 2 , 

computed using the same grid and ^g na i as before. We construct a scatter plot by plotting h^ on one 
and on the other axis, obtaining a 2D plot as the one in Fig.pl As each h u takes values from the 
interval [0,1], with the values being zero (more precisely ~ O(10 )) only for non-resonating sets, the 
plots in Fig. [3] report about qualitative distribution of periodic sets for four computed frequencies. Each 
branch of the scatter plot vectors for a given frequency (defined by a color) represents one periodic chain 
of periodicity integer multiple of the basic period. As the frequencies considered in Fig. [3] are mutually 
prime, the structure of the entire scatter plot approximates the full periodic structure of the phase space 
at e = 0.12 (in the approximation of four frequencies). As opposed to scatter plots of usual time averages 
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Figure 3: The scatter plot for f l {x, y) — sin(77rx + ll7ry) vs. f 2 (x, y) = sin(57rx + 8ny) for e = 0.12 with 
four different frequencies defined by different colors and indicated in the legend. 



discussed in [1 , these scatter plots are less structured as they show the resonating periodic sets only, 
disregarding the rest of the phase space. 

Recall that each point in a scatter plot here has x-coordinate (y-coordinate) given by the value of 
(h^). By clustering the scatter plot points of a given frequency, one can obtain a new correspondence 
between the (available) colors and the initial grid-points in a way analogous to what discussed in [T], in 
order to obtain a better approximation of periodic sets locations in phase space. Again, it is of interest to 
consider as many functions as possible, preferably with both small and large resolution numbers n&m, in 
order to visualize both local and global phase space features. A phase space partition can be constructed 
this way, that visualizes all the periodic sets of a chosen frequency, and therefore giving more insight into 
the phase space structure then a single-function approach discussed earlier. 

5.2 Visualization of the single-frequency periodic partitions 



Consider N linearly independent functions and their corresponding harmonic time averages {h^, . . . h\ 



Their scatter plot will be enclosed in the interval [0, 1]^ with a branching structure similar to Fig.p 



Observe that each scatter plot point is actually a harmonic time average vector h w (xo, yo) corresponding 
to some grid-point (xQ,ya). 

We introduce a norm on the scatter plot-space by considering the Euclidean norm defined on the 
harmonic time average vectors: 



There are clearly more scatter plot points, even with many functions, having the same (or very similar) 
norm ||h w ||. However, note that the correspondence: 



maintains the key properties of the usual harmonic time average: ||h w || is zero for non-resonating periodic 
sets, it is very small for chaotic points and it is O(10 _1 ) order for resonating periodic sets. On the other 
hand, the scalar field ||h w (a;o, yo)|| will yield better visualization of the frequency-^ periodic sets as it "sums 
up" the visualization done by all considered N functions: it approximates the periodic partition of the phase 





(13) 



{xo,Vo) G grid 



s 



space for the given frequency with the approximation-order given by the number of functions involved. It 
is enough for a periodic chain to resonate with at least one considered function, and it will be visualized 
by the ||h w (xo,?/o)||- The downside of this idea is that separate periodic sets belonging to different periodic 
chains might have the same (or very similar) colors assigned: as we are primarily interested in visualizing 
the periodic chains as a whole, this issue will be neglected here. However, by applying phase space zooms 
(similarly to what discussed in [1]), one can analyze the sub-structure of periodic sets at the desired scale. 

We therefore construct a straightforward clustering algorithm of scatter plot points: let a new coloring- 
value for each grid-point (xo,j/o) be given by H u (xo,yo) defined as: 

H u ,{x ,yo) = -~- — , (14) 

m ^(x ,yo)egridllMzo,2/o)|| 

which is the field Hh^xo, yo) || normalized by its biggest phase space value, and therefore with the values 
between and 1. The field H^^xo^yo) approximates the periodic partition for the frequency w. 

In Fig. [4] we show approximations of period-2 and period-5 periodic partitions done with four functions 
of different resolution numbers, and colored (clustered) using the algorithm described above. Harmonic 
time averages were computed for each function separately with the same grid 800 x 800 and Hg na i = 30, 000. 
Observe the improved visualization with respect to the previously shown plots in Fig.[T] By using diverse 




(a) (b) 

Figure 4: Period-2 (a) and period-5 (b) partitions visualized using f 1 = sin(27ra;+37ry), f 2 = sin(37nc+77n/), 
/ 3 = sin(57rx + 87ry) and / 4 = sin(77ra; + liny) for e — 0.12. Log-scale was used for the field. 

resolution numbers we managed to simultaneously visualize periodic chains at all scales (including the 
integer multiple periodicities). Given the initial grid-point resolution, the pictures in Fig.[4ji&b are rather 
exhaustive in terms of the quantity of visualized period-2 and period-5 phase space subsets. However, 
as already mentioned, the price paid amounts to having red color for different unrelated periodic chains 
(essentially for all of them). But on the other hand, note the clear coloring difference between the resonating 
periodic chains (red and dark red), chaotic zone (light green - yellow), and non-resonating periodic chains 
(blue and dark blue). Period-5 partition in Fig. [4)3 even visualizes the secondary chaotic zone between the 
elliptic period-2 points in the middle. 

5.3 Visualization of the ergodic regions 

As observed in the previous Section, harmonic time averages for irrational frequency do not resonate 
with none of the periodic sets, while still resonate weakly with the chaotic zone points (just as any other 
harmonic time average resonates weakly within the chaotic zone). As shown in Fig.|TJl and examined in 
the Section [4], this can be used for visualization of the chaotic zone. Furthermore, from the discussion in 
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the previous paragraph it follows that more functions can be used in order to improve the visualization, in 
a way analogous to the periodic partition approximation. 

In Fig. [5] we show plots obtained by clustering according to the Eq. (14 1 for three functions (of diverse 
resolution numbers). Each harmonic time average was computed separately with the same grid, but for 
a longer transient of fifj na i = 100, 000 and an irrational frequency. Observe the precision in visualization 




of the chaotic zone, specially for the larger e-value in Fig.[5p. We see not only primar (diffusing) and 
secondary chaotic zones (located around the period-2 chain), but also higher-order ones like around the 
period-3 chain (Fig.j5^i). Also, as opposed to previously shown one- function chaotic zone visualization 
(Fig.[l]l), in these plots the coloration of the chaotic zone is more uniform (yellow) and hence more visible. 
This is again due to more functions producing a better approximation for the final field H^^xq, yo). 

In order to show the visualization properties of harmonic time averages in the context of differing 
between the periodic sets and the chaotic regions, in Fig. [6] we show a histogram of H u (xo,yo) values for 
pictures in Figs.[3|k[5j We use the log-scale for in order to have an analogy with the mentioned pictures. 
We indicate three clearly visible peaks common to all distributions: non-resonating periodic set always have 
values ~ O(10 -9 ), chaotic regions always resonate weakly with values ~ O(10~ 5 ), and the periodic sets in 
the case of the appropriate frequency resonate having the values ~ O(10 _1 ). Observe that both rational 
and irrational frequencies respect the mentioned range of -ff^xoj 2/o) values in relation to the dynamical 
behaviors (for non-clustered values the same properties hold but less precisely). This proves that by 
using more functions we improve the quality of visualization in sense of better pronounced peaks. 

Furthermore, a cut-off can be introduced between these peaks in order to visualize only one preferred 
dynamical zone, corresponding to the examined peak. For instance, a full picture of the regular phase 
space portion (all frequencies) could be constructed by cutting out the resonant peaks for many diverse 
rational frequencies/functions and considering corresponding phase space points. Alternatively, the same 
phase space portion could be observed if one considers many functions with irrational frequencies, focusing 
on completely non-resonating part of the phase space. 
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Figure 6: Histograms of values used in Figs. [4]^5] with ui and e values indicated in the legend. Log-scale 
was used for x-axis in analogy with the previous pictures. Three zones with different dynamical behaviors 
and different resonating regimes are indicated. 



6 Extended Standard Map 

In this Section, we show the application of the exposed visualization method on the Extended Standard Map 
proposed in |24] and already studied in [1]. It is a three-dimensional measure-preserving action- action- angle 
map (a particular generalization of the standard map) defined as: 

x' = .t + esin(27rz) + (5sin(27ry) [modi] 

y' — y + esin(27rz) [modi] (15) 



z 



= z + x + esin(27rz) + 5 sin(27ry) [mod 1] 



It has been suggested that this map is ergodic for any non-zero values of perturbations e and 8, backed by 
the observation that no invariant surface persists for any non-zero perturbation [24 . 

We follow the argumentation from |T| and add a further argument to the mentioned claim. We set a 
three-dimensional grid of 20 x 20 x 20 points and evaluate harmonic time averages for two functions with 
the frequencies of u = |, g, g and the transient of ^fi na j — 10 7 . We then cluster the values separately 



for each frequency following the exposed procedure Eq. (14) and plot the histograms for H UJ (x , Vo, z ) 
in Fig. [7] Two different examples of perturbation values are considered: there seems to be only a small 
difference in the profiles that resemble log-normal distributions, while they both tend to zero for all the 
investigated frequencies. This means that at the considered non-zero perturbation values there appears to 
be no persisting periodic sets of any periodicity, which confirms the proposed claim within the given range 
of precision. Moreover, the range of values for the distribution of H u] {xo,yo,Zo) (around 10 -4 ) indicates 
presence of the chaotic dynamics, in the sense of a weak resonance for all the frequencies. We expect that 
for a longer time-evolution the distributions shrink towards zero, confirming the initial ergodic hypothesis. 
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Figure 7: Histograms of H u (xq, yo, zq) obtained by clustering the values of harmonic time averages for 
f 1 = sin(27rx + 3iry + 5ttz) and / 2 = sin(ll7ra; + 6iry + 8irz), for the frequencies ui = 5,3,5 as indicated. 
The case of s = 0.01, J = 0.001 in (a) and e = 0.2, <5 = 0.05 in (b). 



7 Conclusions 

We exposed a new method of frequency analysis of dynamical systems based on harmonic time averages 
3 , following the ergodic theory visualization technique previously proposed in [TJ [5] • Using the known 
properties of standard map, we showed how the method can be implemented in the context of discrete- 
time dynamical systems. Harmonic time averages for various frequencies were computed over the standard 
map's phase space and their precisions estimated. A simple algorithm suggesting a clustering method for 
improving visualization by using more functions was constructed and employed. The periodic partitions of 
different periodicities were correctly visualized, including the corresponding chaotic zones. The distribution 
of values of H u was analyzed, and the peaks corresponding to different types of dynamics were identified. 
The method was applied to the Extended Standard Map, confirming the previous suggestions regarding 
its dynamical properties. 

The most clear improvement of the method lies in the optimization of the clustering algorithm dis- 
cussed in Section [5] A better version should design a parametrization for all the scatter plot branches, 
in order have a continuous correspondence between the scatter plot vectors and the (available) colors. 
This would moreover solve the issue of same-coloring for the independent periodic sets. Furthermore, the 
same discussion from [1] applies in the regard of optimizing the shape/volume of the clustering cells: in 
cases of more complex dynamical systems (higher dimensionality) where scatter plots are not showing the 
branching structure, cell division of scatter plot space might be the only choice. 

In addition, similarly to the case in PQ, relationships between the convergence properties and the 
nature of the underlying orbit is to be investigated. While the convergence slope gives rough estimate of 
the orbit's type, a detailed analysis of the convergence pattern for rational and irrational frequencies might 
yield further insights. 

The method was here exposed using the standard map, as it is a very investigated case of 2D map. 
However, the full applicational extent of the method lies in high-dimensional systems, where 2D phase 
space sections can be considered, in analogy with the discussion in pQ. Furthermore, the method is fully 
applicable to continuous-time dynamical systems as well (ODEs), although the computation of harmonic 
time averages in these cases might be somewhat more demanding. Finally, it would be interesting to see if 
the method is applicable to complex systems like multi-dimensional coupled maps on networks [551 US] . 
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